Examples
In [1]:
Copied!
from bcn.binding_network import binding_network
from bcn.catalysis_network import catalysis_network, binding_and_catalysis
from bcn.binding_network import binding_network
from bcn.catalysis_network import catalysis_network, binding_and_catalysis
In [2]:
Copied!
import numpy as np
import pandas as pd
import time
import sympy as sp
from scipy.integrate import solve_ivp
from scipy.linalg import null_space
import itertools
import cvxpy as cp
import plotly.graph_objects as go
import plotly.express as px
import numpy as np
import pandas as pd
import time
import sympy as sp
from scipy.integrate import solve_ivp
from scipy.linalg import null_space
import itertools
import cvxpy as cp
import plotly.graph_objects as go
import plotly.express as px
Define binding networks¶
There are three reaction networks need to idenfify
- $E + S \leftrightharpoons C$
- $E + 2S \leftrightharpoons C$
- $E + S \leftrightharpoons C_1,\ S + C_1 \leftrightharpoons C_2$
In [3]:
Copied!
def binding_catalysis_network(k_bind, k_cat, bn, total_E, total_const_idx, xcat_in_total_idx, cat_active_in_xbind_idx):
"""
Using Binding and Catalysis Networks (BCN) to smiluate the speed of reaction
inputs:
outputs:
"""
s_mat = np.array([
[-1], # S
[1], # P
])
cn = catalysis_network(s_mat)
kbind = np.array(k_bind) # kes (um)
kcat = np.array(k_cat) # kcat_s2p
bcn = binding_and_catalysis(bn, cn, kbind, kcat, np.array([total_E,]), total_const_idx, xcat_in_total_idx, cat_active_in_xbind_idx)
return bcn
def binding_catalysis_network(k_bind, k_cat, bn, total_E, total_const_idx, xcat_in_total_idx, cat_active_in_xbind_idx):
"""
Using Binding and Catalysis Networks (BCN) to smiluate the speed of reaction
inputs:
outputs:
"""
s_mat = np.array([
[-1], # S
[1], # P
])
cn = catalysis_network(s_mat)
kbind = np.array(k_bind) # kes (um)
kcat = np.array(k_cat) # kcat_s2p
bcn = binding_and_catalysis(bn, cn, kbind, kcat, np.array([total_E,]), total_const_idx, xcat_in_total_idx, cat_active_in_xbind_idx)
return bcn
define binding network: $E + S \leftrightharpoons C$¶
In [4]:
Copied!
# binding part for 1st reaction network
# E + S = C
n_mat=np.array([
[1, 1, 0, -1] # binding reaction E + S <-> C
])
# x_sym = sp.symbols('x_E, x_S, x_P, x_C,')
# t_sym = sp.symbols('q_{E}, q_{S}, t_{P},')
# k_sym = sp.symbols('K_{ES},')
# bn1 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym)
bn1 = binding_network(n_mat, is_atomic=True)
bn1.l_mat
# binding part for 1st reaction network
# E + S = C
n_mat=np.array([
[1, 1, 0, -1] # binding reaction E + S <-> C
])
# x_sym = sp.symbols('x_E, x_S, x_P, x_C,')
# t_sym = sp.symbols('q_{E}, q_{S}, t_{P},')
# k_sym = sp.symbols('K_{ES},')
# bn1 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym)
bn1 = binding_network(n_mat, is_atomic=True)
bn1.l_mat
Out[4]:
array([[ 1., 0., 0., 1.],
[ 0., 1., 0., 1.],
[ 0., 0., 1., -0.]])
define binding network: $E + 2S \leftrightharpoons C$¶
In [5]:
Copied!
# binding part for 2nd reaction network
# E + 2S = C
# binding part
n_mat=np.array([
[1, 2, 0, -1] # binding reaction E + 2S <-> C
])
# x_sym = sp.symbols('x_E, x_S, x_P, x_C,')
# t_sym = sp.symbols('t_{E}, t_{S}, t_{P},')
# k_sym = sp.symbols('K_{ES},')
# bn2 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym)
bn2 = binding_network(n_mat, is_atomic=True)
bn2.l_mat
# binding part for 2nd reaction network
# E + 2S = C
# binding part
n_mat=np.array([
[1, 2, 0, -1] # binding reaction E + 2S <-> C
])
# x_sym = sp.symbols('x_E, x_S, x_P, x_C,')
# t_sym = sp.symbols('t_{E}, t_{S}, t_{P},')
# k_sym = sp.symbols('K_{ES},')
# bn2 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym)
bn2 = binding_network(n_mat, is_atomic=True)
bn2.l_mat
Out[5]:
array([[ 1., 0., 0., 1.],
[ 0., 1., 0., 2.],
[ 0., 0., 1., -0.]])
define binding network: $E + S \leftrightharpoons C_1,\quad S + C_1 \leftrightharpoons C_2$¶
In [6]:
Copied!
# binding part for 3rd reaction network
# E + S = C1, S + C1 = C2
n_mat = np.array([
[1, 1, 0, -1, 0],
[0, 1, 0, 1, -1],
])
# x_sym = sp.symbols('x_E, x_S, x_P, x_{C1}, x_{C2},')
# t_sym = sp.symbols('q_{E}, q_{S}, q_{P},')
# k_sym = sp.symbols('K_{ES}, K_{SC1},')
# bn3 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym)
bn3 = binding_network(n_mat, is_atomic=True)
bn3.l_mat
# binding part for 3rd reaction network
# E + S = C1, S + C1 = C2
n_mat = np.array([
[1, 1, 0, -1, 0],
[0, 1, 0, 1, -1],
])
# x_sym = sp.symbols('x_E, x_S, x_P, x_{C1}, x_{C2},')
# t_sym = sp.symbols('q_{E}, q_{S}, q_{P},')
# k_sym = sp.symbols('K_{ES}, K_{SC1},')
# bn3 = binding_network(n_mat, is_atomic=True, x_sym=x_sym, t_sym=t_sym, k_sym=k_sym)
bn3 = binding_network(n_mat, is_atomic=True)
bn3.l_mat
Out[6]:
array([[ 1., 0., 0., 1., 1.],
[ 0., 1., 0., 1., 2.],
[ 0., 0., 1., -0., -0.]])
Define catalysis networks¶
define catalysis network: $S \rightarrow P$¶
In [7]:
Copied!
s_mat = np.array([
[-1], # S
[1], # P
])
cn = catalysis_network(s_mat)
s_mat = np.array([
[-1], # S
[1], # P
])
cn = catalysis_network(s_mat)
Define binding and catalysis networks¶
In [8]:
Copied!
total_E = 5 # total enzyme
total_S = 5 + 1e-6
total_P = 1e-6
total_E = 5 # total enzyme
total_S = 5 + 1e-6
total_P = 1e-6
define BCN1¶
bn1 and cn
In [9]:
Copied!
k_bind_1 = np.array([1]) # kbind = k- / k+
k_cat_1 = np.array([1]) # kcat_s2p
# bcn
total_const_idx_1 = np.array([0,]) # qe is first in totals (qe, qs, qp)
xcat_in_total_idx_1 = np.array([1, 2]) # qs is 2nd element in totals, tp is not contained.
cat_active_in_xbind_idx_1 = np.array([3]) # c_es
k_bind_1 = np.array([1]) # kbind = k- / k+
k_cat_1 = np.array([1]) # kcat_s2p
# bcn
total_const_idx_1 = np.array([0,]) # qe is first in totals (qe, qs, qp)
xcat_in_total_idx_1 = np.array([1, 2]) # qs is 2nd element in totals, tp is not contained.
cat_active_in_xbind_idx_1 = np.array([3]) # c_es
In [10]:
Copied!
bcn1 = binding_and_catalysis(
bn1,
cn,
k_bind_1,
k_cat_1,
np.array([total_E,]),
total_const_idx_1,
xcat_in_total_idx_1,
cat_active_in_xbind_idx_1
)
bcn1 = binding_and_catalysis(
bn1,
cn,
k_bind_1,
k_cat_1,
np.array([total_E,]),
total_const_idx_1,
xcat_in_total_idx_1,
cat_active_in_xbind_idx_1
)
define BCN2¶
bn2 and cn
In [11]:
Copied!
# binding & catalysis network for 2nd reaction network
# E + 2S = C
k_bind_2 = np.array([1]) # kbind = k- / k+
k_cat_2 = np.array([1]) # kcat_s2p
# bcn
total_const_idx_2 = np.array([0,]) # qe is first in totals (qe, qs, qp)
xcat_in_total_idx_2 = np.array([1, 2]) # qs is 2nd element in totals, tp is not contained.
cat_active_in_xbind_idx_2 = np.array([3]) # c_es
# binding & catalysis network for 2nd reaction network
# E + 2S = C
k_bind_2 = np.array([1]) # kbind = k- / k+
k_cat_2 = np.array([1]) # kcat_s2p
# bcn
total_const_idx_2 = np.array([0,]) # qe is first in totals (qe, qs, qp)
xcat_in_total_idx_2 = np.array([1, 2]) # qs is 2nd element in totals, tp is not contained.
cat_active_in_xbind_idx_2 = np.array([3]) # c_es
In [12]:
Copied!
bcn2 = binding_and_catalysis(
bn2,
cn,
k_bind_2,
k_cat_2,
np.array([total_E,]),
total_const_idx_2,
xcat_in_total_idx_2,
cat_active_in_xbind_idx_2
)
bcn2 = binding_and_catalysis(
bn2,
cn,
k_bind_2,
k_cat_2,
np.array([total_E,]),
total_const_idx_2,
xcat_in_total_idx_2,
cat_active_in_xbind_idx_2
)
define BCN3¶
bn3 and cn
In [13]:
Copied!
# binding & catalysis network for 2nd reaction network
# E + S = C1, S + C1 = C2
k_bind_3 = np.array([1, 1]) # kbind = k- / k+
k_cat_3 = np.array([0.1, 0.1]) # kcat_s2p
# k_cat_3 = np.array([1, 1])
# k_cat_3 = np.array([0.5, 0.5])
# bcn
total_const_idx_3 = np.array([0,]) # qe is first in totals (qe, qs, qp)
xcat_in_total_idx_3 = np.array([1, 2]) # qs is 2nd element in totals, tp is not contained.
cat_active_in_xbind_idx_3 = np.array([3, 4]) # c_es
# binding & catalysis network for 2nd reaction network
# E + S = C1, S + C1 = C2
k_bind_3 = np.array([1, 1]) # kbind = k- / k+
k_cat_3 = np.array([0.1, 0.1]) # kcat_s2p
# k_cat_3 = np.array([1, 1])
# k_cat_3 = np.array([0.5, 0.5])
# bcn
total_const_idx_3 = np.array([0,]) # qe is first in totals (qe, qs, qp)
xcat_in_total_idx_3 = np.array([1, 2]) # qs is 2nd element in totals, tp is not contained.
cat_active_in_xbind_idx_3 = np.array([3, 4]) # c_es
In [14]:
Copied!
bcn3 = binding_and_catalysis(
bn3,
cn,
k_bind_3,
k_cat_3,
np.array([total_E,]),
total_const_idx_3,
xcat_in_total_idx_3,
cat_active_in_xbind_idx_3
)
bcn3 = binding_and_catalysis(
bn3,
cn,
k_bind_3,
k_cat_3,
np.array([total_E,]),
total_const_idx_3,
xcat_in_total_idx_3,
cat_active_in_xbind_idx_3
)
Simulating the dynamic process¶
In [15]:
Copied!
def get_trajs(bn, bnc, k_bind, total_S, total_P, start_time=0, end_time=20, timestamps=20):
logxcat_init = np.log10(np.array([total_S - total_P, total_P])) # S, P
a_mat = bn.l_mat
bnc.kbind = k_bind
bnc.logkbind = np.log10(k_bind)
trajs = {}
time, logxcat_traj, logxbind_traj, logder_traj = bnc.get_traj(logxcat_init, start_time, end_time, timestamps, a_mat=a_mat)
trajs['x_cat'] = 10 ** logxcat_traj
trajs['x_bind'] = 10 ** logxbind_traj
trajs['logder'] = logder_traj
trajs['time'] = time
return trajs
def get_trajs(bn, bnc, k_bind, total_S, total_P, start_time=0, end_time=20, timestamps=20):
logxcat_init = np.log10(np.array([total_S - total_P, total_P])) # S, P
a_mat = bn.l_mat
bnc.kbind = k_bind
bnc.logkbind = np.log10(k_bind)
trajs = {}
time, logxcat_traj, logxbind_traj, logder_traj = bnc.get_traj(logxcat_init, start_time, end_time, timestamps, a_mat=a_mat)
trajs['x_cat'] = 10 ** logxcat_traj
trajs['x_bind'] = 10 ** logxbind_traj
trajs['logder'] = logder_traj
trajs['time'] = time
return trajs
In [16]:
Copied!
idx_P_cn = 1
trajs_1 = get_trajs(bn1, bcn1, k_bind_1, total_S, total_P)
trajs_2 = get_trajs(bn2, bcn2, k_bind_2, total_S, total_P)
trajs_3 = get_trajs(bn3, bcn3, k_bind_3, total_S, total_P)
fig = go.Figure()
fig.add_trace(go.Scatter(x=trajs_1['time'], y=trajs_1['x_cat'][:, idx_P_cn].tolist(),
mode='lines', name='bn1',
))
fig.add_trace(go.Scatter(x=trajs_2['time'], y=trajs_2['x_cat'][:, idx_P_cn].tolist(),
mode='lines', name='bn2',
))
fig.add_trace(go.Scatter(x=trajs_3['time'], y=trajs_3['x_cat'][:, idx_P_cn].tolist(),
mode='lines', name='bn3',
))
fig.update_layout(
font=dict(
size=27
),
template='simple_white',
width=800, height=400,
)
fig.show(renderer='notebook')
idx_P_cn = 1
trajs_1 = get_trajs(bn1, bcn1, k_bind_1, total_S, total_P)
trajs_2 = get_trajs(bn2, bcn2, k_bind_2, total_S, total_P)
trajs_3 = get_trajs(bn3, bcn3, k_bind_3, total_S, total_P)
fig = go.Figure()
fig.add_trace(go.Scatter(x=trajs_1['time'], y=trajs_1['x_cat'][:, idx_P_cn].tolist(),
mode='lines', name='bn1',
))
fig.add_trace(go.Scatter(x=trajs_2['time'], y=trajs_2['x_cat'][:, idx_P_cn].tolist(),
mode='lines', name='bn2',
))
fig.add_trace(go.Scatter(x=trajs_3['time'], y=trajs_3['x_cat'][:, idx_P_cn].tolist(),
mode='lines', name='bn3',
))
fig.update_layout(
font=dict(
size=27
),
template='simple_white',
width=800, height=400,
)
fig.show(renderer='notebook')
In [17]:
Copied!
trajs_1['x_cat'][:, idx_P_cn].shape
trajs_1['x_cat'][:, idx_P_cn].shape
Out[17]:
(20,)
Visulize three ROP space¶
In [18]:
Copied!
logmin = -6
logmax = 6
chart='x'
def sample(bn, b_vec, kcat_list=None, logmin=-6, logmax=6, npoints=100, chart='x', logvar=None):
# randomly sample values for x
logvar = np.random.uniform(logmin, logmax, (npoints, bn.dim_n)) if logvar is None else logvar
# calculate log derivatives
logders, logx = bn.logder_num(logvar, chart=chart)
if isinstance(b_vec, list):
dots_in_rop = [bn.logder_activity_num(b, logx, logders) * k for b, k in zip(b_vec, kcat_list)]
dots_in_rop = sum(dots_in_rop)
else:
dots_in_rop = bn.logder_activity_num(b_vec, logx, logders)
return dots_in_rop, 10 ** logx
logmin = -6
logmax = 6
chart='x'
def sample(bn, b_vec, kcat_list=None, logmin=-6, logmax=6, npoints=100, chart='x', logvar=None):
# randomly sample values for x
logvar = np.random.uniform(logmin, logmax, (npoints, bn.dim_n)) if logvar is None else logvar
# calculate log derivatives
logders, logx = bn.logder_num(logvar, chart=chart)
if isinstance(b_vec, list):
dots_in_rop = [bn.logder_activity_num(b, logx, logders) * k for b, k in zip(b_vec, kcat_list)]
dots_in_rop = sum(dots_in_rop)
else:
dots_in_rop = bn.logder_activity_num(b_vec, logx, logders)
return dots_in_rop, 10 ** logx
In [19]:
Copied!
bcn1_ROP, bcn1_x_points = sample(bn1, b_vec=np.array([0, 0, 0, 1]))
bcn2_ROP, bcn2_x_points = sample(bn2, b_vec=np.array([0, 0, 0, 1]))
bcn3_ROP, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, k_cat_3[0], k_cat_3[1]]))
bcn3_ROP_1, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, 1, 0]))
bcn3_ROP_2, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, 0, 1]))
bcn1_ROP, bcn1_x_points = sample(bn1, b_vec=np.array([0, 0, 0, 1]))
bcn2_ROP, bcn2_x_points = sample(bn2, b_vec=np.array([0, 0, 0, 1]))
bcn3_ROP, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, k_cat_3[0], k_cat_3[1]]))
bcn3_ROP_1, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, 1, 0]))
bcn3_ROP_2, bcn3_x_points = sample(bn3, b_vec=np.array([0, 0, 0, 0, 1]))
In [20]:
Copied!
E_index = 0
S_index = 1
P_index = 2
E_index = 0
S_index = 1
P_index = 2
In [21]:
Copied!
ROP = go.Figure()
ROP.add_trace(go.Scatter(x=bcn1_ROP[:, S_index], y=bcn1_ROP[:, E_index],
mode='markers', name=r'bn1',
))
ROP.add_trace(go.Scatter(x=bcn2_ROP[:, S_index], y=bcn2_ROP[:, E_index],
mode='markers', name=r'bn2',
))
ROP.add_trace(go.Scatter(x=bcn3_ROP[:, S_index], y=bcn3_ROP[:, E_index],
mode='markers', name=r'bn3 combination',
))
ROP.add_trace(go.Scatter(x=bcn3_ROP_1[:, S_index], y=bcn3_ROP_1[:, E_index],
mode='markers', name=r'bn3 C1',
))
ROP.add_trace(go.Scatter(x=bcn3_ROP_2[:, S_index], y=bcn3_ROP_2[:, E_index],
mode='markers', name=r'bn3 C2',
))
ROP.update_layout(
font=dict(
size=18
),
template='simple_white',
width=560, height=420,
xaxis_title=r'$\Large{\frac{\partial \log x_C}{\partial \log q_S}}$',
yaxis_title=r'$\Large{\frac{\partial \log x_C}{\partial \log q_E}}$',
xaxis=dict(range=[-1, 3]),
yaxis=dict(range=[-1, 3], scaleanchor = "x", scaleratio = 1),
)
ROP.show(renderer='notebook')
ROP = go.Figure()
ROP.add_trace(go.Scatter(x=bcn1_ROP[:, S_index], y=bcn1_ROP[:, E_index],
mode='markers', name=r'bn1',
))
ROP.add_trace(go.Scatter(x=bcn2_ROP[:, S_index], y=bcn2_ROP[:, E_index],
mode='markers', name=r'bn2',
))
ROP.add_trace(go.Scatter(x=bcn3_ROP[:, S_index], y=bcn3_ROP[:, E_index],
mode='markers', name=r'bn3 combination',
))
ROP.add_trace(go.Scatter(x=bcn3_ROP_1[:, S_index], y=bcn3_ROP_1[:, E_index],
mode='markers', name=r'bn3 C1',
))
ROP.add_trace(go.Scatter(x=bcn3_ROP_2[:, S_index], y=bcn3_ROP_2[:, E_index],
mode='markers', name=r'bn3 C2',
))
ROP.update_layout(
font=dict(
size=18
),
template='simple_white',
width=560, height=420,
xaxis_title=r'$\Large{\frac{\partial \log x_C}{\partial \log q_S}}$',
yaxis_title=r'$\Large{\frac{\partial \log x_C}{\partial \log q_E}}$',
xaxis=dict(range=[-1, 3]),
yaxis=dict(range=[-1, 3], scaleanchor = "x", scaleratio = 1),
)
ROP.show(renderer='notebook')